Glass Transition of the Monodisperse Gaussian Core Model 
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We numerically investigate the dynamical properties of the one-component Gaussian core model 
in supercooled states. We find that nucleation is increasingly suppressed with increasing density. 
The system concomitantly exhibits glassy, slow dynamics characterized by the two-step stretched 
exponential relaxation of the density correlation and a drastic increase of the relaxation time. We 
also find a weaker violation of the Stokes-Einstein relation and a smaller non-Gaussian parameter 
than in typical model glass formers, implying weaker dynamic heterogeneities. Additionally, the 
agreement of the simulation data with the prediction of mode-coupling theory is exceptionally good, 
indicating that the nature of the slow dynamics of this ultra-soft particle fluid is mean-field-like. 
This fact may be understood as a consequence of the long-range nature of the interaction. 

PACS numbers: 64.70.pv, 64.70.dg,81.05.Kf, 61.43.Fs 



The nature of the glass transition is surrounded by con- 
troversy. Several scenarios have been proposed to explain 
the drastic slowing down of dynamics of supercooled flu- 
ids near the glass transition point [ll-Q . Numerical simu- 
lation of simple model fluids is an ideal tool to test these 
competing scenarios. However, the typical model flu- 
ids studied so far, such as Lennard- Jones, soft-core, and 
hard-sphere mixtures, have short-ranged, strong repul- 
sive interactions in common, which dictate their thermo- 
dynamic, structural, and dynamical properties and ren- 
der the results of these models qualitatively similar 
A new class of model glass formers is desirable to di- 
versify our pictures and perspectives on the glass tran- 
sition within the limited accessible time windows of the 
simulations. Recently, ultra-soft particle fluids have at- 
tracted particular attention in soft-materials science Q. 
They are systems composed of spherical particles inter- 
acting with bounded and weak repulsions and are a good 
model for various soft materials, such as star-polymers 
and dendrimers. The absence of the hard-core-like repul- 
sion makes the thermodynamic and dynamic behaviors of 
this class of systems extremely rich compared with stan- 
dard molecular systems. Their phase diagrams exhibit 
exotic and counterintuitive properties, including a stable 
fluid phase at high temperatures for arbitrary densities, 
re-melting of solids at higher densities, and complex crys- 
talline phases at low temperatures The dynamics of 
the ultra-soft particles fluids also exhibits rich and non- 
trivial behaviors 

In this Letter, we consider the simplest version of ultra- 
soft particles, i.e., the Gaussian core model (GCM) fluid 
originally introduced by Stillinger [lj| . The GCM inter- 
action is given by 
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where e and a characterize the energy and length scales, 
respectively. The GCM is an ideal model to study 
glassy dynamics because its thermodynamic phase di- 
agram is relatively simple. Other ultra-soft particles, 



such as Hertzian spheres and star-polymers, exhibit com- 
plex crystalline phases, which may affect the dynamics in 
the supercooled state 0, |Tl| . We numerically study the 
monodisperse GCM in three dimension and show that nu- 
cleation is suppressed at very high densities and that the 
system exhibits canonical glassy behavior in the super- 
cooled state. The quantitative agreement of the dynam- 
ical properties with the theoretical predictions is better 
than that of all previously investigated model glass for- 
mers. 

The thermodynamic and dynamic properties of the 
GCM have recently been vigorously investigated 0, H, 
Il2l - [l3 |. Most previous studies, however, focused on the 
density regime pa 3 < 1, where the phase diagram ex- 
hibits reentrant melting. At these densities, the monodis- 
perse GCM easily nucleates to form crystals as it crosses 
the phase boundary. In this Letter, we investigate dy- 
namics of the GCM near the fluid-crystal phase boundary 
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FIG. 1: GCM phase diagram (empty squares). Results of 
Prestipino et al. [T2] (filled circles) are also shown. The 
dashed line is a fit by logT m oc -p 2/3 0. The melting 
and freezing lines are indistinguishable at this scale. Crosses 
denote the state points where the MD simulations are per- 
formed. 
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FIG. 2: Time dependence of the orientational order pa- 
rameter Q 6 . (a) p = 0.5, T — 2.5 x 10~ 3 , (b) p = 1.0, 
T = 2.5 x 10 -4 , and (c) p = 1.5, T = 2.4 x 10~ 5 . £ is 
scaled by r Q , which is evaluated from non-crystallizing sam- 
ples. A = (T m —T)/Tm defines the distance from the melting 
temperature. 



at the unprecedentedly high densities of per 3 > 1. 

The thermodynamic properties of the system at high 
densities are carefully characterized using a Monte Carlo 
(MC) simulation. We identify the fluid-crystal phase 
boundary using a thermodynamic integral calculation 
combined with the particle-insertion method and the 
Frenkel-Ladd technique (I2T Uij , as shown in Figure [T] 
Stillinger showed that the ground state of the GCM at 
per 3 > 0.18 is the bec crystal and argued that the melting 
temperature T m obeys log T m oc — p 2 / 3 at high densities 
based on the duality relation with the hard-sphere sys- 
tem [lj| • We find that T m follows this scaling at per 3 > 1 
and confirm that the crystalline structure is indeed bec 
at these densities. The details of the thermodynamic 
properties will be discussed in a forthcoming paper fl6| . 

Dynamics of the system is investigated using a molec- 
ular dynamics (MD) simulation in the NVT ensemble 
with a Nose thermostat in the cubic cell with a peri- 
odic boundary condition. A time-reversible integrator, 
similar to the velocity- Verlet method, is used with a po- 
tential cut-off at r — 5a [jlj . In the following, we take 
a, e/fce and e^m/e) 1 / 2 as the length, temperature, and 
time units, respectively. We focus on the four densities, 
p = 0.5, 1.0, 1.5, and 2.0 (the melting temperatures are 
T m = 4.4 x 10~ 3 , 5.0 x 10~ 4 , 5.8 x 10~ 5 , and 7.2 x lO" 6 , 
respectively) and perform the MD simulations for various 
temperatures below T m , indicated by crosses in Figure [T] 
For each state point, five independent runs are performed 
to improve the statistics. The system size is fixed at 
N = 3456. The simulations for N = 2000 and 9826 con- 
firm that the finite-size effect is negligible. Starting from 
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FIG. 3: F s (k,t) for several state points, (a) p = 0.5 and 
T x 10 3 = 7, 4, 3, 2.6, (b) p = 1 and T x 10 4 = 7, 4, 3, 2.6, 
2.5, (c) p = 1.5 and T x 10 5 = 7, 4, 3, 2.6, 2.4, 2.3, and (d) 
p = 2 and T x 10 6 = 10, 5, 4, 3.4, 3.2, 3, 2.93. The dashed 
lines in (a)-(c) denote the lowest-temperature data for which 
at least one of trajectories crystallizes. The dash-dotted lines 
in (d) are the solutions of the MCT equation. 



the initial configurations generated at high temperatures, 
we perform the simulations for longer than 50r Q , where 
r a is the alpha-relaxation time obtained from the inter- 
mediate scattering function (see below). The nuclcation 
of the system into the crystalline state is monitored by 
the orientational order parameter Qq (l7j . Qq is known 
to be 0.5 for the bec crystal and zero for the disordered 
or fluid state 17J . Figure [2] shows the time evolution of 
Qq for several trajectories at three state points. To com- 
pare the different states on equal footing, t is scaled by 
r a , which is a good measure of the equilibration time. At 
p = 0.5 and T = 2.5 x 10 -3 , for which the distance from 
the phase boundary defined by A = (T m —T)/T m is about 
0.43, all five trajectories crystallize near t ~ 20tq, (Figure 
[2] (a)). At p = 1.0 and T = 2.5 x 10" 4 , despite deeper 
supercooling (A = 0.5) than that in Figure Ufa), nucle- 
ation takes place at much longer times of roughly 40r a 
(Figure [2] (b)). For higher densities, like p = 1.5 (Figure 
[5] (c)), all sampled trajectories remain in the fluid state 
although the system is more supercooled (A 0.6). In 
addition to Q§, we also monitor the potential energy of 
the system, which discontinuously decreases as the sys- 
tem crystallizes. We have checked that its time evolu- 
tion is synchronized with Q§. From these observations, 
we conclude that the nucleation of the GCM is strongly 
suppressed at very high densities. 

Next, we focus on the slow dynamics of the samples 
that do not crystallize. We evaluate the self-part of the 
intermediate scattering function, F s (k, t), after equilibra- 
tion runs of about 30r a . Figure [3] shows the results for 
F s (k,t) for the wavevector k near the first peak of the 
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FIG. 4: (a) The radial distribution function for p = 1.5 and 
2.0 at the lowest temperatures. Ri„t represents the "inter- 
action range" (see text), (b) The MCT power-law fit of r Q 
where T c is a fitting parameter, (c) r^-dependence of the SE 
relation and (d) the peak value of the NGP a m ax- Data for 
the small and large particles of the KA model are also plotted. 
The MCT results are shown in dashed lines. 



static structure factor S(k) at various temperatures for 
four densities. Dashed lines at the lowest temperatures in 
Figure[3](a)-(c) represent the data for the state points at 
which at least one of the five trajectories crystallizes. At 
the lowest density p = 0.5, crystallization always takes 
place before the slow dynamics sets in. At higher densi- 
ties, however, the system clearly exhibits two-step relax- 
ation, while the radial distribution functions g{r) remain 
liquid-like, as shown in Figure|4](a). The sudden appear- 
ance of a plateau in F s (k,t) is the hallmark of the slow 
dynamics near the glass transition. The alpha-relaxation 
time, r a , defined by F s (k, r Q ) = e^ 1 , drastically increases 
as the temperature decreases. In particular, glassy relax- 
ation is observed up to the lowest accessible temperature 
with no indication of crystallization for the densest sys- 
tem p = 2.0. We also calculate the mean-square displace- 
ment (AR 2 (t)) and observe the typical glassy behavior 
characterized by a plateau followed by diffusive behavior 
in the alpha-relaxation regime and the drastic decrease 
of the self-diffusion coefficient D. 

We make a more detailed characterization of the slow 
dynamics and compare the results with the prediction 
of mode-coupling theory (MCT). MCT has successfully 
described many dynamical properties of moderately su- 
percooled fluid using S(k) as the sole input Q. Though 
still contentious, MCT is believed to be a dy namic mean- 
field theory of the glass transition [l|, Il8l [lif . It pre- 
dicts the relaxation behaviors of correlation functions 



such as F s (k,t) semi-quantitatively and the power-law 
increase of r a ~ \T — T^l -7 , where T c is the temper- 
ature at which MCT predicts the spurious nonergodic 
transition. Other properties that MCT successfully pre- 
dicts include the time-temperature superposition (TTS) 
in the alpha-relaxation regime, the fe-dependence of the 
plateau height of the intermediate scattering function, 
and dynamic scaling in the plateau regime [3j. On the 
other hand, MCT fails to capture dynamics below T c , 
where the activation processes over the complex energy 
landscape dominate. Another failing of MCT is that the 
T c 's obtained by the fitting of simulation data systemat- 
ically deviate from those evaluated from the theory fioj ]. 
Furthermore, due to the mean-field nature of the theory, 
MCT lacks an explanation of the violation of the Stokes- 
Einstcin (SE) relation and growth of non-Gaussian pa- 
rameters (NGP) (2l|. We solve the MCT equation for 
the GCM using S(k) obtained from simulation and com- 
pared the solution with the simulation data. We mainly 
focus on the data for p = 1.5 and p = 2.0, for which 
the plateau of the two-step relaxation of F s (k,t) is well 
developed. Our results confirm that the F s (k,t) sim- 
ulation data obeys TTS in the alpha-relaxation regime 
and can be fitted by a stretched exponential function 
e -(t/T„)' 3 w j|;h t ne exponent j3 rs 0.7, which agrees with 
the MCT results. The fc-dependence of the plateau height 
of F s (k, t) agrees with MCT as well. We also find that the 
temperature dependence of r Q follows the MCT power 
law, r a oc \T — j , c ( sim )|-7 j as shown in Figure|4] (b), where 
r a is plotted using the short-time relaxation time to de- 
fined by F s (k,to) = 0.95 as a time unit. The result for 
the binary Lennard- Jones system (KA model) [22[ is also 
shown. Wc fit the data using 7 sa 2.7 (for both p = 1.5 
and 2.0), which is obtained from the MCT solution, and 
leaving Tc S as the fitting parameter. Note that the 



value of 7 is comparable to that of the KA model [22 [. 
These results demonstrate that the GCM shares many 
properties with other model glass formers. As shown be- 
low, however, the GCM is distinct in several respects. 
First, the agreement of MCT's F s (k,t) with the simula- 
tion data is very good. The dash-dotted lines in Figure [3] 
(d) are the MCT solution fitted using e = 1-T/T C as the 
sole parameter (aside from the time unit). This agree- 
ment is striking, given that for other model fluids e (and 
sometimes the wavevectors as well) needs to be adjusted 
at each temperature to obtain a reasonable fit [2(| [23| (an 
exception is the four-dimensional system (24|). Second, 
the parameters T , j snn ' 1 used to fit r Q in Figure |4] (b) are 



unprccedentedly close to the theoretical values T c 



(theory) 



We find that 



(sim) 



2.02 x 10~ 5 and 2.66 x 10" 6 for 



p = 1.5 and 2.0, respectively, whereas their theoretical 
counterparts are j , c ( thcory ) = 2.66 x 10~ 5 (p = 1.5) and 
3.20 x 10~ 6 (p = 2.0). The discrepancies between the 
simulation and theory are 32% and 20% for p = 1.5 and 
2.0, respectively. For other glass formers, Tj slm ) i s known 



4 



to differ considerably from T, 



(theory) 



for example, T, 



(theory),. 



For the KA model, 



0.92 compared to T c 



(sim) 



0.44. 



The discrepancy is more than 100% [20(. The KA model 
is still a high-temperature fluid and F a (k,t) 



(theory) 



at T t 

decays exponentially without a hint of two-step relax- 
ation. On the contrary, the GCM at j , c ( theory ) ij es Jeep in 
the region where the plateau of F s (k,t) is well developed 
(see FigurclHJd)). The third and most noticeable point is 
the suppression of the SE violation and the NGP. In Fig- 
ure |4] (c) and (d), Dr a normalized by high-temperature 
values (-DTa)ref and the peak value of the NGP a max are 
shown as a function of r a /t (instead of T to compare 
the different systems on equal footing). The NGP is de- 
fined by a (t) = 3(Ai? 4 (i)) /5{AR 2 (t)} -1. The variations 
of both Dr a and a max for the GCM are much weaker 
than those of the KA model [22| ■ Similar suppression of 
the SE violation was observed in four-dimensional sys- 
tems 24, 2f|. Furthermore, a max is smaller for p = 2.0 
than for 1.5. Because the SE violation and the growth 
of the NGP are thought to be the consequences of the 
underlying dynamic heterogeneities near the glass tran- 
sition point [2l| , our results imply that the dynamic het- 
erogeneities are weaker in the GCM and, thus, the nature 
of glassy dynamics of the GCM is more mean-field-likc 
than those of other systems [28| . 

The mean-field nature of the GCM may be attributed 
to the long-range nature of the interaction potential at 
high densities and extremely low temperatures, where 
many particles interact with each other. In Figure |4] (a) , 
"interaction range" Rint defined by v(r = Ri n t) = k^T 
is indicated by dashed lines. Rint reaches the second 
and third coordination shells, which means that many 
particles enter in the range of Rint at these high densi- 
ties. This is in stark contrast with ordinary fluid systems 
with strong repulsive interactions for which the interac- 
tion range is on the order of a or the distance between 
neighboring particles. A more detailed analysis support- 
ing this hypothesis is reported elsewhere [16j . 

An explanation of the drastic decrease of nucleation 
rates of the GCM at high densities is still lacking. It 
is tempting to speculate that this phenomenon is inti- 
mately related to the mean-field character of the GCM. 
In the context of classical nucleation theory, the time 
scale of nucleation r„ is proportional to that of transla- 
tional diffusion td ~ 1/D. Recently, Tanaka has argued 
that nucleation should always intervene before the dy- 
namic arrest takes place if the SE relation is violated 
because the decoupling of the translational motion of a 
single particle and structural relaxation leaves r„ insensi- 
tive to temperature, while the bulk dynamics drastically 
slow down [2^|. This scenario has been recently exam- 
ined numerically [22]] ■ The opposite may take place for 
the GCM; that is, the weaker SE violation may lead to 
the concomitant increase of r„ and T a , ultimately sup- 
pressing nucleation. 



In conclusion, we demonstrated that the GCM is an 
unexpectedly simple and novel glass former. The rich 
dynamics of the GCM and ultra-soft particle systems in 
general may answer some important unanswered ques- 
tions regarding the glass transition and nucleation. 
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sions. A. I. is supported by the JSPS and K. M. by 
KAKENHI No. 2154016 and Priority Areas "Soft Mat- 
ter Physics" . 
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